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data,  one  of  the  standard  regularization  methods(see  Poggio  and  Torre,  1984) 
leads  to  cubic  spline  interpolation  before  differentiation.  We  show  that 
in  the  case  of  regularly-spaced  data  this  solution  corresponds  to  a  convolution 
filter —  to  be  applied  to  the  signal  before  differentiation —  which  is  a 
cubic  spline.  In  the  case  of  non-exact  data,  which  is  the  most  interesting 
situation,  we  use  another  regularization  method  that  leads  to  a  different 
variational  principle.  We  prove(l)  that  this  variational  principle  leads 
to  a  convolutioal  filter  for  the  problem  of  one  dimensional  edge  detection, 

(2)  that  the  form  of  this  filter  is  very  similar  to  the  gaussian  filter, 

and  (3)  that  the  regularizing  parameter  A  in  the  variational  principle  effectively 

controls  the  scale  of  the  filter. 

Finally,  we  outline  several  issues  arising  from  our  solution  to 
the  edge  detection  problem:  (1)  the  use  of  methods  from  regularizing  theories 
for  finding  the  optimal  value  of  the  regularizing  parameter  A  ;  (2)  the  connection 
between  these  methods  and  the  scale-space  method  for  edge  detection;  (3) 
the  relationship  between  our  edge  detector  and  other  detectors,  especially 
the  Marr/  Hildreth  edge  detector;  (4)  the  extension  of  our  one-dimensional 
solutionto  two-dimensional  edge  detection;  and  (5)  the  extension  of  our  method 
to  deal  with  differentiation  of  surface  data(  though  the  physical  constraint 
underlying  the  form  of  the  regularizer  is  not  valid  in  general  for  depth 
data);  this  issue  is  connected  to  the  problem  of  interpolating  and  approximating 
depth  data. 
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Abstract.  -We  consider  edge  detection  as  the  problem  of  measuring  and  localizing  changes 
of  light  intensity  in  the  image.  /As  discuss^^ by  Torre  and  Poggio  (1984),  edge  detection, 
when  cefined  in  this  way,  is  an  llTposed  problem  in  the  sense  of  Hadamafd.  —  — — 

Using  standard  regularization  theory,  we  regularize  the  problem  with  a  stabilizing 
functioial  that  is  a  specific  form  of  a  Tikhonov  stabilizer,  following  Reinsch  (1967)  and 
Schoerberg  (1964).  The  regularized  solution  that  arises  is  then  the  solution  to  a  variational 
principle.  In  the  case  of  exact  data,  one  of  the  standard  regularization  methods  (see  Poggio 
and  Torre,  1984)  leads  to  cubic  spline  interpolation  before  differentiation.  We  show  that  in 
the  case  of  regularly-spaced  data  this  solution  corresponds  to  a  convolution  filler--to  be 
applied  to  the  signal  before  differentiation— which  is  a  cubic  spline.  In  the  case  of  non-exact 
data,  which  is  the  most  interesting  situation,  we  use  another  regularization  method  that  leads 
to  a  different  variational  principle.  We  prove  (1)  that  this  variational  principle  loads  to  a 
convolution  filter  for  the  problem  of  one  dimensional  edge  detection,  (2)  that  the  form  of 
this  filter  is  very  similar  to  the  gaussian  filter,  and  (3)  that  the  regularizing  parameter  X  in 
the  variational  principle  effectively  controls  the  scale  of  the  filter. 

Finally,  we  outline  several  issues  arising  from  our  solution  to  the  edge  detection 
problem:  (1)  the  use  of  methods  from  regularizing  theories  for  finding  the  optimal  value 
of  the  regularizing  parameter  X;  (2)  the  connection  between  these  methods  and  the  scale- 
space  method  for  edge  detection;  (3)  the  relationship  between  our  edge  detector  and 
other  detectors,  especially  the  Marr/Hildreth  edge  detector;  (4)  the  extension  of  our  one¬ 
dimensional  solution  to  two-dimensional  edge  detection;  and  (5)  the  extension  of  our  method 
to  deal  with  differentiation  of  surface  data  (though  the  physical  constraint  underlyi  ig  the 
form  of  the  regularizer  is  not  valid  in  general  for  depth  data);  this  issue  is  connected  to  the 
problem  of  interpolating  and  approximating  depth  data. 
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1.  Introduction 


^  Edge  detection  does  not  have  a  precisely  defined  goal.  The  word  edge  itself,  which  refers 
to  physical  properties  of  objects,  is  somewhat  of  a  misnomer.  Several  years  of  experience 
have  shown  that  the  ideal  goal  of  detecting  and  locating  physical  edges  in  the  surfaces 
being  imaged  is  very  difficult  and  still  out  of  reach  ^{fi^  a  review  see  Brady,  198^  Edge 
detection  has  come  to  be  defined  as  the  first  step  in  this  goal  of  detecting  physical  changes 
such  as  object  boundaries— the  operation  of  detecting  and  locating  changes  in  intensity  in 
the  image.  Other  processes  which  operate  on  these  measurements  of  intensity  changes  will 
then  group  boundaries  and  label  and  characterize  them  in  terms  of  the  properties  of  the 
3-D  surfaces. 

Intended  in  this  narrow  sense,  edge  detection— this  first  step  in  processing  the 
image — is  mainly  the  process  that  measures,  detects  and  localizes  changes  of  intensity. 
Derivatives  must  be  estimated  correctly  to  label  the  critical  points  in  the  image  intensity 
array,  characterize  their  local  properties  (are  they  minima  or  maxima  or  saddle  points?) 
and  thi  s  relate  them  to  the  underlying  physical  process  (are  they  shadow  edges  or  depth 
discont  nuities?)./cAs_a_consequence,  several  different  derivatives  of  the  image,  possibly  at 
different  scales,  may  have  to  be  esfifftated.^ - - — _ _ 

In  this  sense,  Torre  and  Poggio  (1984)  considered  edge  detection  as  a  protJem  of 
numerical  differentiation  of  images.  The  problem  is  not  straightforward,  and  attempts  over 
many  years  have  proven  its  difficulties.  Considered  as  a  problem  of  numerical  differentiation, 
edge  detection  turns  out  to  be  an  ill-posed  problem.  As  explained  by  Poggio  and  Torre 
(1984),  mathematically  ill-posed  problems  are  problems  where  the  solution  either  does  not 
exist  or  is  not  unique  or  does  not  depend  continuously  on  the  data. 

Numerical  differentiation  is  a  (mildly)  ill-posed  problem  because  its  solution  does  not 
depend  continuously  on  the  data.  It  is  therefore  natural  to  try  to  solve  this  problem  by  using 
regularization  techniques  developed  in  recent  years  for  dealing  with  mathematically  ill-posed 
problems.  The  problem  can  be  regularized  by  the  use  of  a  wide  class  of  filters  (Torre  and 
Poggio,  1984,  section  2.4;  see  also  Duda  and  Hart,  1973).  In  the  following  section  we 
consider  two  specific  regularizing  operators,  that  in  some  sense  are  very  natural. 

2.  Regularizing  Edge  Detection 

To  regularize  an  ill-posed  problem  and  make  it  well-posed,  one  has  to  introduce  generic 
constraints  on  the  problem.  In  this  way,  one  attempts  to  force  the  solution  to  lie  in  a  subspace 
of  the  solution  space,  where  it  is  well  defined.  The  basic  idea  of  regularization  techniques 
is  to  restrict  the  space  of  acceptable  solutions  by  choosing  the  function  that  minim  zes  an 
appropriate  functional.  Poggio  and  Torre  consider  in  particular  standard  regularization 
theory  based  on  quadratic  variational  principles.  They  list  three  main  techniques  for 
regularizing  the  ill-posed  problem  of  finding  z  from  the  data  y  such  that  Az  =  y  They 
involve  the  choice  of  norms  H  H  (usually  quadratic)  and  of  a  stabiiizing  functional  l|7’*||- 
The  choice  is  dictated  by  mathematical  considerations,  and,  most  importantly,  by  a  physical 
analysis  of  the  generic  constraints  on  the  problem.  Three  main  methods  can  then  be  sipplied 
(see  Bertero,  1982): 

(1)  Among  ?  that  satisfy  ||P«1|  <  C,  where  C  is  a  constant,  find  z  that  minimizes 

||>i*-v||,  (1) 


’  A  very  similar  problem  arises  in  the  characterization  of  surface  properties— in  particular  their 
differential  properties— from  depth  data. 


Poggio 


Voorhees,  Yuille 


Figure  1  Cubic  Interpolating  Spline  Filter  (L4  function  of  Schoenberg) 


(2)  Among  z  that  satisfy  \\Az  -  j/j|  <  C,  find  z  that  minimizes 

\\P4.  (2) 

(3)  Find  2  that  minimizes 

||>l2-,,||*  +  X||P2||»,  (3) 

where  X  is  a  regularization  parameter. 

The  first  method  consists  of  finding  the  function  2  that  satisfies  the  constraint  \\Pz\\  <  C  and 
best  approximates  the  data.  The  second  method  computes  the  function  2  that  is  sufficiently 
close  to  the  data  (C  depends  on  the  estimated  errors  and  is  zero  if  the  data  are  noiseless) 
and  is  most  “regular”.  In  the  third  method,  the  regularization  parameter  X  controls  the 
compromise  between  the  degree  of  regularization  of  the  solution  and  its  closeness  to  the 
data. 

In  the  case  of  edge  detection  considered  as  numerical  differentiation,  we  want  an 
approximation  /  to  the  intensity  data  yi  at  sample  points  x,  that  is  well  behaved  under 
differentiation.  Thus  we  consider  an  operator  A  which  samples  the  function  /  on  the  lattice 
such  that  A/I*.  =  /I*,  for  t  = 

The  problem  is  then  to  find  a  suitable  norm  and  a  suitable  stabilizing  functional  ||P/|l. 
It  is  natural  to  chose  for  P  the  simplest  form  of  Tikhonov's  stabilizing  functionals  (Tikhonov 
and  Arsenin,  1976)  with  P  =  ^  and  the  usual  L*  norm.  For  fc  =  2  this  choice  corresponds 

to  a  constraint  of  smoothness  on  the  approximated  intensity  profile  2,  with  ||F’/||  = 

Its  physical  justification  is  that  the  noiseless  image  has  to  be  smooth  in  the  sense  that  its 
derivatives  must  be  bounded  because  the  image  is  band-limited  by  the  optics.  Band-limited 
functions  have  bounded  derivatives  because  /'  <  fiM,  where  M  =  sup  F{w),  n  is  the 
cut-off  frequency,  and  F(u)  is  the  Fourier  transform  of  /(x).  Physically,  the  constraint  of 
smoothness  allows  us  to  effectively  eliminate  the  noise  that  creeps  in,  after,  or  during  the 
sampling  and  transduction  process  and  makes  the  operation  of  differentiation  unstable.  We 
stress  that  this  is  not  the  only  stabilizing  functional  possible  for  this  problem,  although  it  is 
probably  the  simplest  one. 

2.1.  The  second  regularization  method  and  interpolating  cubic  splines 

With  this  choice  of  P,  the  second  regularization  method  is:  among  /  such  that  /(x.)  =  y, 
find  /  that  minimize  //"’  dx.  A  theorem  by  Schoenberg  (see  Greville,  1969)  shows  that 
the  solution  to  this  problem  is  a  cubic  spline  interpolating  between  the  data  points.  The 
following  result  is  a  reformulation  of  results  due  to  Schoenberg  (1946,  1964); 

Theorem  1:  The  cubic  spline  function  f  interpolating  evenly-spaced  data  points 
and  minimizing  J  J"‘  dx  can  be  obtained  by  convolving  the  data  points  with  a 
cubic  spline  filler  which  corresponds  to  the  L4  function  of  Schoenberg. 

A  plot  of  Lt  is  given  in  Figure  1.  Note  that  the  size  of  the  filter  is  fixed  with  respect  to 
the  sampling  lattice.  The  filter  is  0  at  every  pixel  but  the  central  one,  where  it  is  1 .  Thus  L4 
is  an  interpolating  filter  that  does  not  perform  any  smoothing. 
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Figure  2  Filter  R  derived  by  regularization  principles. 


2.2.  T:ie  third  regularization  method,  approximating  splines  and  the  edge  de '.action 
filter 

The  third  regularization  method  leads  to  the  following  problem:  Find  the  /  that  minimizes 

iV'  ~  /(^‘))*  +  ^  (4) 

i 

where  is  the  regularization  parameter  which  can  be  found  as  described  later.  This  problem 
was  considered  originally  by  Reinsch  (1967)  in  the  case  of  numerical  differentiation  and  by 
Schoenberg  (1964)  for  the  problem  of  graduation.  Both  Schoenberg  and  Reinsch  gave 
the  soljtion  in  terms  of  approximating  cubic  splines.  In  addition,  we  prove  that  far  most 
practicil  purposes,  the  approximating  spline  function  can  be  obtained  by  convohing  the 
data  pciint  y,  with  the  cubic  spline  convolution  filter  R  shown  in  Figure  2  (see  also  Appendix 
1).  We  then  have  the  following  theorem: 

Theorem  2:  The  solution  to  equation  (4) — the  reqularized  solution  to  the  problem 
of  numerical  differentiation— in  the  case  of  inexact  data,  can  be  obtained  by 
convolving  the  data  with  a  convolution  filter  which  is  (a)  a  cubic  spline,  and  ’b) 
very  similar  to  a  gaussian. 

Although  this  result  is  especially  significant  in  the  context  of  edge  detection  where 
the  search  for  an  optimal  filter  and  its  justification  has  been  a  longstanding  preoccupation, 
it  is  somewhat  surprising  that  this  result  does  not  seem  to  have  been  widely  appieciated 
in  the  numerical  analysis  literature.  The  exact  assumptions  under  which  Theonjm  2  is 
valid  are  discussed  in  Appendix  1.  First,  the  data  must  be  given  on  a  regular  grij  (as  is 
the  case  for  an  image)  Second,  the  image  data  must  either  go  to  zero  at  infinity  or  be 
periodic.  Under  these  conditions,  the  filtering  operation  is  space-invariant  and  linear  (the 
Euler-Lagrange  equations  corresponding  to  the  quadratic  variational  problem  are  linear). 
Thus  the  approximating  spline  can  be  obtained  by  a  convolution  operation.  Note  ‘hat  the 
result  that  the  regularizing  operator  corresponding  to  a  quadratic  variational  principle  is  a 
convolution  filter— for  data  on  a  regular  grid  and  toroidal  boundary  conditions— is  /alid,  in 
general,  beyond  the  case  of  numerical  differentiation. 

Theorem  3:  Quadratic,  Tikhonov  type  regularization  principles  are  equivalent  to 
convolving  the  data  with  a  generalized  spline  filter,  if  the  data  are  given  on  an 
regular  lattice  and  the  boundary  conditions  are  appropriate. 

A  generally  interesting  question  is  the  physical  correctness  of  the  regularized  folution. 
In  the  c  ase  considered  in  this  paper  the  answer  is  simple  and  not  very  insightful:  a  necessary 
and  sufficient  condition  for  the  regularized  solution  to  be  correct  is  that  the  true  intensity 
distribution  is  a  polynomial  of  order  less  than  4  between  sampling  points.  This  property  can 
be  derived  directly  from  equation  (3)  of  Appendix  1.2. 

3.  Regularization  parameter  and  comparison  with  the  gaussiari  filter 

Figure  2  shows  the  filter  R  obtained  by  solving  the  variational  principle  equatioi  (1)  in 
Appendix  11.  Its  shape  and  size  depends  on  the  regularizing  parameter  X.  Figure  3(a) 
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Figure  3  Filter  R  and  first  derivative  for  different  values  of  regularizing  parameter  X.  (a) 
X  affects  size  of  (the  first  derivative  of)  R  but  (b)  does  not  appear  to  affect  shape  of 
R.  (Amplitudes  are  normalized  in  (a);  both  amplitudes  and  widths  are  scaled  linearly  and 
independently  in  (b)  for  comparison.) 


Figure  4  Comparison  of  one-dimensional  regularizing  filter  R  with  Gaussian:  (a)  zeroth, 
(b)  first,  and  (c)  second  derivatives. 

shows  the  first  derivative  of  the  filter  for  different  values  of  X.  The  continuous  version  of 
the  filter,  derived  in  Appendix  2,  is  practically  indistinguishable  from  the  discrete  filter,  as 
shown  by  numerical  comparisons. 

It  is  rather  intuitive  that  the  smoothing  parameter  X  controls  the  effective  size  of  the 
filter.  From  our  numerical  work,  it  seems  that  X  does  not  significantly  affect  the  shape  of 
the  filter,  but  only  its  size,  as  shown  in  Figure  3(b).  Changing  X  amounts  to  scaling  the  size 
of  the  filter  up  or  down.  If  X  is  small,  smoothness  is  unimportant,  and  the  fiiter  will  tend  to 
be  an  interpolating  filter  and  therefore  be  similar  to  a  ^  function.  On  the  other  hand,  with 
a  very  large  X,  the  main  weight  is  on  smoothness,  and  the  filter  will  tend  to  be  very  large. 
The  continuous  form  of  the  filter  suggests  that  the  role  of  X  is  indeed  equivalent  to  the  role 
of  a  for  the  gaussian  (X  ~  aV  as  shown  in  Appendix  4). 

The  regularization  filter  derived  here  appears  to  be  quite  similar  to  the  Gaussian 
distribution.  Graphs  of  the  filter  R,  its  first  and  second  derivatives  are  shown  with  those 
of  the  Gaussian  in  Figure  4.  Marr  and  Hildreth  (1980;  Hildreth,  1980)  have  argued  that 
the  Gaussian  is  an  optimal  smoothing  filter  for  detection  due  to  its  localization  properties 
in  both  the  spatial  and  frequency  domains.  The  fact  that  the  Gaussian  is  quite  similar  to 
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the  optimal  filter  derived  here  using  regularization  principles  provides  further  mathf'matical 
justification  for  the  use  of  a  Gaussian-like  filter  in  edge  detection. 

If  the  boundary  conditions  are  not  periodic  or  natural,  then  the  derivation  of  l^einsch 
(see  Appendix  1.1)  provides  the  correct  Green  function  for  those  boundary  corditions. 
In  numerical  experiments  we  have  found  that  the  Green  function  obtained  in  this  way  is 
space-invariant  for  all  points  but  those  close  to  the  boundaries. 

4.  Discussion 

Severe'  questions  and  extensions  suggest  themselves  in  a  natural  way.  Here  we  li  3t  some 
of  them. 

4.1.  F  nding  the  Optimal  X 

Regularization  theory  can  give  the  optimal  value  of  X  if  the  errors  on  the  smoothing  criteria 
and  th'j  error  of  the  approximations  are  known  in  advance.  If  the  integral  of  f"‘  is  less 
than  £,  and  the  sum  of  (k,  -  /(xi))^  is  less  than  t,  then  X  =  (/£’  (see  Bertero,  1981; 
Tikhon  3v,  1963;  Tikhonov  and  Arsenin,  1977).  Normally,  however,  errors  on  the  data  or 
on  the  smoothness  conditions  are  not  known  in  advance.  Regularization  theory  provides 
several  methods  for  finding  the  optimal  smoothing  parameter  X  under  this  circumstance. 
We  want  to  indicate  here  two  main  methods:  (1)  Tikhonov's  method,  for  convolut.on  type 
problems,  as  is  the  case  here,  and  (2)  the  cross  validation  method  and  the  generalized 
cross-validation  method  (Wahba,  1980). 

We  plan  to  evaluate  these  methods  for  finding  the  optimal  X  for  edge  detection  and 
to  test  them  on  real  images.  An  interesting  issue  that  we  are  also  planning  to  evpiore  is 
the  following:  The  basic  idea  of  the  generalized  cross-validation  method  is  to  check  the 
goodness  of  approximation  for  each  value  of  X.  In  order  to  do  that,  one  computes  the 
approximation  by  using  not  all  but  only  some  of  the  data  points.  Thus,  the  data  po  nls  that 
are  not  used  for  computing  the  approximations  serve  as  the  control  points  for  the  goodness 
of  the  approximation.  If  one  computes  the  goodness  of  the  fit  at  different  X,  one  can  then 
choose  the  optimal  X.  This  idea  has  obvious  connections  with  the  use  of  fingerprints  (Yuille 
and  Poggio,  1983)  for  finding  the  natural  scale  of  the  filter  (Witkin,  1983);  this  point  is 
discussed  next. 

4.2.  Optimal  X  and  natural  scale 

The  size  of  the  filter  with  which  to  perform  edge  detection  has  always  been  an  uniesolved 
issue  in  computer  vision.  Our  approach  makes  it  clear  that  one  expects,  indeed,  an  optimal 
size  of  the  filter  associated  with  the  optimal  value  of  the  smoothing  parameter  X.  In  more 
recent  years,  several  scales  of  filtering  have  been  used,  partly  as  a  way  around  this  problem. 
Rosenfeld  and  Kak  (1982),  Marr  and  Poggio  (1977),  Marr  and  Hildreth  (1980;  Hildreti,  1980) 
have  used  several  sizes  of  filters  in  order  to  perform  edge  detection. 

More  recently,  Witkin  (1983)  has  suggested  the  use  of  scale-space  filtering,  essentially 
filtering  across  a  continuum  of  scales,  as  the  method  by  which  to  choose  the  optim.il  scale. 
Witkin  suggested  some  heuristics  for  picking  the  natural  scale  of  filtering.  We  belisve  that 
cross  validation  type  methods  may  make  more  rigorous  the  idea  of  selecting  an  optimal 
filterinc  scale  by  selecting  the  optimal  X  value.  We  are  planning  to  use  fingerprints — the 
zero-crossings  across  scales  of  the  Laplacian  of  the  convolved  image— to  find  the  optimal 
X, 

For  scale  space  filtering  to  be  maximally  effective,  the  shape  of  the  filter  should  be 
a  gaussian  (Yuille  and  Poggio,  1983a,  1983b;  Witkin,  1983).  The  effect  of  changing  X  is 
essentially  equivalent  to  changing  the  size  of  the  filter,  and  furthermore,  the  underly  ng  filter 
is  very  similar  to  a  gaussian.  Therefore,  increasing  the  value  of  X  is  equivalent  to  filtering 
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Figure  5  Difference  of  two-dimensional  Gaussians  as  approximations  to  (a)  Laplacian  of 
gaussian,  (b)  Laplacian  of  two-dimensional  regularizing  filter  Ri.  The  ratio  between  the 
scales  (a)  of  the  two  Gaussians  is 


the  signal  with  a  larger  and  larger  gaussian  as  required  in  scale  space  filtering.  Because 
of  the  many  nice  properties  of  gaussian  convolution,  this  also  suggests  an  efficient  way  of 
obtaining  approximations  at  increasing  X. 

4.3.  Relation  with  other  edge  detectors 

Because  of  the  close  similarity  of  our  cubic  spline  filter  to  a  gaussian,  the  edge  detector 
that  we  derive  in  this  paper  is  very  similar  to  edge  detectors  proposed  previously.  Marr 
and  Poggio  (1977)  proposed  the  difference  of  two  gaussians  as  an  approximation  to  the 
second  derivative  of  a  gaussian.  Marr  and  Hildreth  (1980;  Hildreth,  1980)  have  shown  that 
the  second  derivative  of  a  gaussian  is,  indeed,  very  close  to  the  difference  of  gaussians, 
J.  Canoy’s  filter  (Canny,  1983)  is  very  close  to  the  derivative  of  a  gaussian,  and  Haralick’s 
cubic  polynomial  interpolant  (Haralick,  1982)  is  again  similar  to  Canny's  filter. 

Our  derivation  justifies  the  use  of  a  gaussian  or  a  filter  very  close  to  a  gaussian  as 
the  best  filter  for  edge  detection.  Regularization  theory  yields  derivative-of-gaussian  filters 
as  the  optimal  filter  in  a  simpler,  more  general,  and,  we  believe,  more  rigorous  way,  than 
previous  derivations.  In  particular,  our  result  makes  clear  that  the  quasi-gaussian  filter 
regularizes  the  ill-posed  problem  of  numerical  differentiation.  The  regularizing  constraint 
here  is  that  the  norm  of  the  derivatives  in  the  noise-free  image  is  small. 

It  is  interesting  that  we  derive  a  filter  very  similar  to  Canny's,  based  on  simpler  and 
more  general  principles  that  are  not  restricted  to  the  optimal  detection  of  step  edges.  It  is 
also  interesting  to  note  that  the  second  derivative  of  the  regularization  filter,  like  the  second 
derivative  of  the  Gaussian,  can  be  approximated  by  a  difference  of  Gaussians  (although 
not  as  well).  While  the  second  derivative  of  the  Gaussian  is  best  approximated  by  a  space 
constant  ratio  (the  ratio  of  scales  of  the  two  Gaussians)  ~t  ~  l.r>  (Marr  and  Hildreth  1980; 
Hildreth,  1980),  increasing  the  ratio  to  7  =s  I  results  in  a  function  which  better  fits  the  main 
(excitatory)  lobe  of  the  regularizing  filter,  as  shown  in  Figure  5. 

4.4.  Extension  to  two  dimensions 

Our  approach  can  be  extended  to  two  dimensions  in  several  ways.  The  most  straightforward 
method  involves  the  use  of  directional  derivatives.  First  order  directional  ilerivatives  are 
taken  along  several  directions.  Fach  one  of  them  is  one-dimensional  and  can  be  performed 
according  to  our  one-dimensional  edge  detector  scheme.  Do()ending  on  the  specific  goal, 
one  may  then  choose  the  direction  that  gives  the  maximum  value  of  the  derivative.  This 
corresponds  to  the  non-ma.ximum  supression  scheme  used  by  Canny  (1003).  It  is  also 
equivalent  to  taking  the  second  directional  derivative  along  the  gradient  and  looking  for  its 
zeros  (Torre  and  Poggio,  1984). 
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Figure  6  Cross  section  of  (a)  two-dimensional  regularization  filter  U-i  and  (b)  its  Laplacian. 


A  second  method  requires  directly  formulating  the  regularization  principle  in  two 
dimensions.  Instead  of  equation  (2),  one  would  then  have  the  problem  of  minimizing 

+  ^  /  /  (/'  ^2/  (>'») 

The  main  problem  is  the  choice  of  the  operator  /'.  If  we  consider  a  Tikhonov  stabilizer  (see 
Poggio  and  Torre,  1984),  then  a  choice  for  /’  that  is  smooth  enough  to  allow  the  use  of 
second  derivatives  of  the  regularized  image  is  /'  V^V.  This  choice  is  used  in  Appendix 
3  to  derive  the  filter  for  the  two-dimensional  continuous  case.  The  filter  is  shown  in  Figure 
6.  The  choice  of  the  derivative  to  be  used  on  the  filter  is  a  separate,  important  issue  that 
we  do  not  address  in  this  paper.  Torre  and  Poggio  (1984)  discuss  the  properties  of  several 
two-dimensional  differential  operators,  including  the  second  directional  derivative  aicng  the 
gradient. 

If  /’  is  chosen  to  be  the  quadratic  variation  or  the  square  Laplacian,  the  resulting 
approximations,  knov;n  as  as  thin  plate  splines  (Wahba,  1980;  Terzopoulos,  1984a),  are 
not  smooth  enough  for  finding  zeros  of  second  derivatives  of  a  function^,  as  implied  by 
Terzopoulos  (1984b).  It  may  also  be  interesting  to  explore  the  filters  resulting  rom  a 
non-linear  functional  /’.  In  the  case  of  a  non-linear  one  cannot  use,  in  general,  the 
standard  results  about  uniqueness  and  other  properties  of  the  solution  that  are  available  for 
the  quadratic  case  of  Tikhonov’s  stabilizers,  because  the  functional  is  no  longer  convex. 

Clearly,  formulations  of  tliis  type  are  also  relevant  for  the  problem  of  surface  interpola¬ 
tion  and  approximation  in  the  sense  of  Gri.mson  (1982)  and  Terzopoulos  (1984a).  In  the 
case  of  sparse  data,  which  they  considered,  the  variational  principle  does  not  lead  to  a 
convolution  filter,  although  it  does  lead  to  a  standard  Green  function.  On  a  regular  grid  it 
leads  to  a  convolution  filter  similar  to  the  gaussian.  As  a  practical  implication,  evenly-spaced 
suiface  data  (for  example,  laser  range  data)  may  be  interpolated  or  approximated  effe  ctively 
by  gaussian  convolution.  Hence,  tasks  wfiich  involve  differentiating  surface  data,  such 
as  con>puting  lines  of  curvature  (Brady,  Ponce,  Yuille  and  Asada,  1984),  could  use  the 
simpler  convolution  method  to  smooth  the  data.  Since  Reinsch's  method  (see  Appendix  1.1) 
can  deil  with  f)oundary  conditions  different  from  periodic  ones,  the  corresponding  Green 
(unction  can  be  used  to  prevent  smoothing  across  depth  discontinuities. 

Tl  e  results  of  applying  the  Laplacian  of  the  two-dimensional  regularization  filtar  and 
the  l.ar.  lacian  of  a  gai.issian  to  an  image  are  shown  in  Figure  7.  As  expected,  due  to  the 
similaridy  of  the  two  filters,  both  edge  detection  operators  yield  similar  results. 


'Wo  are  indcr.'li'd  tii  Dotnotri  Tcr/opoulos  for  this  remark 
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Figure  7  Comparison  of  Laplacian  of  two-dimensional  regularization  filter  and  Laplacian  of 
Gaussian  as  Edge  Detectors;  (a)  image  /,  (b)  zero-crossings  of  (c)  zero- crossings 

of 
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Appendix  1:  The  regularization  filter  for  discrete  data  in  one  dimension 

Below  we  present  two  methods  for  proving  that  the  third  regularization  method  described 
in  the  text  yields  a  convolution  filter  in  the  case  of  evenly-spaced  discrete  data  n  one 
dimension  with  appropriate  boundary  conditions.  Each  derivation  includes  a  method  for 
computing  the  filter,  which  has  the  form  of  a  cubic  spline. 


1.1.  The  convolution  version  of  Reinsch’s  optimal  filter 


Reinsch  (1967)  considered  the  general  problem  of  interpolating  and  smoothing  a  sequence 
of  (not  necessarily  evenly-spaced)  data  points  with  uncertainties.  Given  a  sequence  of 
points  (y.)  with  uncertainties  iAj/,),  representing  values  of  a  function  at  points  (x,),  where 
I  =  1, 2  . . ri  and  I,  <  I i  <  <  the  problem  of  finding  a  smooth  function  which  aasses 

near  each  point  can  be  formulated  as  finding  the  function  /(i)  that  minimizes  the 

functional  (to  be  compared  with  equation  4  in  the  text) 


/(j.)  -  Vx 

^Vi 


(/"(x))*dx. 


(1) 


Using  calculus  of  variations,  a  cubic  spline  function  is  shown  to  be  the  optimal  function 
satisfying  (1), 

!{x]  =  a,  ^  b,(x  -  c,(i- Xif +  di{x-Xif,  (2) 

for  I  =-  1, 2, . ,  n  1 .  and  formulas  are  presented  for  calculating  from  the  data  (x,)  and  (y,) 
the  coefficients  (a,),  (6,),  (c,),  and  (d,). 

W?  specialize  Reinsch's  derivation  as  follows:  First,  the  uncertainty  associated  with 
each  data  point  is  assumed  to  be  the  same  so  that  all  ^,  =  1.  Equation  (1)  above  thus 
reduces  to  equation  (4)  in  the  text.  Second,  the  data  are  assumed  to  be  uniformly  spaced, 
so  that  X,.. I  -  j,  =  /i  for  t  =  i,2,...,n  -  1,  Finally,  unlike  Reinsch’s  derivation,  we  assume 
periodic  boundary  conditions,  /.e.,  that  the  solution  /  is  a  periodic  function  over  !?,  with 
/(i)  =  /(x  ±  tn)  for  integer  k. 

With  these  specializations  we  show  that  each  set  of  coefficients  can  be  calculated  by 
multipymg  the  data  points  y  by  a  constant  coefficient  matrix: 


ft  =  c  —  Cy,  i  =  Dy, 

where  A.  Ih  C,  and  D  are  circulant  matrices  representing  the  filter  and  its  derivatives. 
Hence,  the  data  can  be  optimally  smoothed  (and  differentiated)  by  a  convolution  operation. 

Using  calculus  of  variations,  one  obtains  from  (1)  these  conditions  for  the  optimal 
function  f(x): 

/(3^.)+ - =  0  (3a) 

/'(x.)+  -  /'(x.)-  =  0  (36) 

f"(x,U  -  /"(x.)-  =  0  (3c) 

r[x,U-r(x,).  =  -i{f{x,)-y,),  (3rf) 


(using  =  1)  and 


/""(x)  =  0,  X,  <  X  <  xi+i, 


(4) 


for  «  —  1 _ ,n  (with  xi  —  i„+i),  where 


/<*’(x,)±  =  lim  ±  c). 

«— 0 
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Equations  (JM)  show  that  /(i)  is  a  cubic  spline,  having  the  form  of  (2);  hence, 

/(l,)+  =  Oi 
/'(*<)+  =  *. 

/''(i,)+  =  2c, 

nx.)+  =  8d. 

and,  since  all  i,  -  *,_i  =  A, 


f{xi)-  —  fli-i  +  bi-ih  +  c,-_ifc*  +  di-ih* 

/'(xj)-  =  6,-1  +  2c,_iA  +  3(i<_iA* 

/"(x.)_  =  2c._,  +  6(f._iA 
nx.)_  =  6rf.„. 

Now  substitution  into  (3)  yields 

Oj  —  o,_i  -  —  c<_iA*  —  di-ih*  =  0 

6,-  —  Ai— 1  —  2c, -—lA  —  3<i,'_iA*  =  0 
2c,-  —  2e,-_|  —  6<f,'_iA  =  0 

idi  -  Mi-i  =  -  j^(/(x,)  -  y,). 

Equations  (7)  can  be  manipulated  to  yield 

6, A  =  (a,+j  -  o,)  -  c,A®  -  dih* 
a<+i  —  2aj  +  a,-_i  =  j(c,-_i  +  4c,-  +  c,-+i)A* 
dih*  —  |(c,-+i  -  c,)A* 

(c,-_i  -  2c,  +  Ci+i)A®  —  ^{yi  -  o,)A* 

Using  the  notation  that; 


/  =  the  n  X  n  identity  matrix  [«,-,*]  where 


t  „  =  P’  'O  =  fc,  j  = 

\0,  otherwise 

iV  =  the  n  X  n  “next”  matrix  [n,,*]  where 

„  /l.  if  j  =  A- Imodn,  y  =  1 . n 

\o,  otherwise 
P  =  the  n  X  n  "previous”  matrix  = 


we  define  the  matrices 


For  example,  if  n  =  4, 


1  0 
0  ] 
0  0 
.0  0 


7’=|P+|/+jAr 
Q  =  P-2I+N. 


0  1  0  0\ 
0  0  1  o| 
0  0  0  1 
>1  0  0  0/ 


N  = 


0  0  0  1\ 
1  0  0  o| 
0  1  0  ol 
.0  0  10/ 


4/3 

1/3 

0 

1/3 


1/3  0  1/3 

4/3  1/3  0 

1/3  4/3  1/3 
0  1/3  4/3, 


\ 

(-2 

1 

0 

1  \ 

1  ^  = 

1 

-2 

1 

0 

0 

1 

-2 

1 

1 

^  1 

0 

1 

-2/ 

(5a) 

(56) 

(5c) 

(5rf) 


{6a) 

(86) 

(6c) 

M 

(7a) 

(76) 

(7c) 

(7d) 


(8a) 

(86) 

(8c) 

m 
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Defining  the  vectors  a  =  (a,)’’,  b  =  (b.h)'’’,  c  =  [cih^f,  d  =  and  g  =  (y,)^,  for 

i  —  1, . .  .,n,  equations  (8)  can  be  express^  as 

b  —  {N  -  I)a-c-d 
Q&=  Ti 
d^\{N-I)i 


These  can  be  simplified  to 

a  =  Ay,  6  =  iiy,  c  =  (7y,  d  = 

(10) 

where 

(11a) 

A  =  I-%QC 

(116) 

D  =  j(N-/)C 

(11c) 

B  =  (JV~f)A-C-D. 

(lid) 

Since  I,  N,  and  P  are  circulant  matrices,  and  since  the  set  of  circulant  matrices  is 
closed  under  matrix  addition,  multiplication,  and  inversion,  the  resulting  coefficient  matrices 
A,  B,  C,  and  D  are  also  circulant,  representing  filters  R(x),  R\x),  R"(x),  and  R"'(x).  The  fact 
that  the  filters  are  derivatives  of  each  other  follow  from  (5);  hence,  since  differentiation  is  a 
linear  operator,  the  filter  R",  for  example,  can  be  used  to  both  optimally  smooth  and  twice 
differentiate  the  data. 

1 .2  Another  derivation  of  the  regularization  filter  for  discrete  data  in  one  dimension 

In  this  derivation  we  show  that  the  third  regularization  method,  described  by  equation  (4) 
of  the  text,  yields  a  convolution  filter  which  is  a  cubic  spline.  Standard  results  from  the 
calculus  of  variations  guarantee  that  our  solution  has  continuous  second  derivatives. 

Again,  the  problem  is  to  minimize 

^  f{f'\^)fdx  +  -  Vi)*-  (1) 

t 

We  find  the  minimum  by  sending  /(i)  /(x)  +  6f(x)  and  setting  the  first  variation  of  (1)  to 

zero, 

X  f  f""(x)  6f{x)dx  +  ^{/{xO  -  Vi)  6f[xi)  =  0.  (2) 

« 

This  yields  the  Euler-Lagrange  equation 

X/""(x)  +  /(x)  ^  S{x  -  I.)  =  ^  Vi6{x  - 1,).  (3) 

i  i 

So  far,  we  have  deliberately  not  specified  boundary  conditions.  For  infinite,  or  toroidal, 
boundary  conditions,  the  function  /(x)  can  be  determined  in  terms  of  the  (y.)  by  a  convolution 
filter  if  and  only  If  the  data  points  (li)  are  evenly  spaced.  This  is  because  the  system  is 
then  translation-invariant^.  We  will  show  this  explicitly  and  give  a  method  for  constructing 
the  convolution  filter. 

The  function  in  (1)  is  convex,  and  hence  has  a  unique  minimum,  so  there  is  a  unique 
solution  for  J(x)  in  (3).  Thus,  we  only  need  to  see  whether  a  convolution  can  solve  (3).  We 


(9a) 

(96) 

(9c) 

(9(0 


^Boundary  conditions  other  than  infinite  or  toroidal  will  destroy  translation  invariance. 


Poggio 


Voorhees.  Yuitle 


try  a  solution 

f{x)  =  «(x)  •  y{x)  =  j  R{x-  (4) 

where  denotes  convolution,  K(i)  is  a  filter,  and  y(i)  =  -  *•)•  We  substitute  (4) 

into  (3)  and  obtain 

^  -  *.)  +  Sm  -  ijX*  -  *.)  -  “  *•)  =  ®  (5) 

I  «•  j  • 

We  compare  coefficients  of  yi  and  obtain 

\R""(x  -  I,)  +  ^  R{xj  -  Xi)S(x  —  Xj)  -  S(x  -  Xi)  =  0  (8) 

} 

If  the  (xi)  are  not  evenly-spaced,  then  these  equations  are  inconsistent,  and  no  convolution 
filter  exists.  If  the  (z,)  are  evenly  spaced,  then  the  set  of  equations  in  (6)  reduces  to  a  single 
equation 

\R""{x)  +  Y.  -  *>)  -  «(*)  =  0  (7) 

J 

The  solutions  to  (7)  correspond  to  cubic  splines  “stitched"  together  at  the  points  {z{}. 
Let  Ri(x)  denote  the  solution  in  the  range  x,  <x  <  z,+,.  We  write 

Ri{x)  =  a,z*  -t-  /9,x*  7,x  +  Si  (8) 

The  splines  are  stitched  together  so  that  R{x),  R'(x),  and  R"(x)  are  continuous  at  the  points 
{i,}.  From  (7),  we  see  that  R'"  has  a  discontinuity  of  -iR(xi)  at  i,.  It  is  straightforward  to 
find  the  relations  between  the  «,(x)  and  Ri-i(x)  In  terms  of  the  parameters  a,,/?,,  7,  and  Si. 
This  gives 

On-H  =  an  - 
^n+l  - 

ln-¥l  —  In  2)^”^ 

6n+l  =  Sn  +  ^■gj“7Z(nfc) 

where  h  is  the  spacing  between  the  lattice  points,  h  =  xi+i-  x,,  and 

R{nh)  —  a„(n/»)*  +  ^„(nh)*  +  7„(n/i)  +  (10) 
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Appendix  2:  The  regularization  filter  for  continuous  data  in  one  dimension 
We  want  to  find  the  function  /(x)  which  minimizes  the  quadratic  functional 

/(/(*)  -  +  >'  /  (/"(*))*«**•  (1) 

where  y(x)  are  the  data,  given  on  the  continuum.  We  will  show  that  /(x)  can  be  expressed 
as  a  convolution  of  the  data  y(x)  with  a  suitable  filter  R(x). 

We  obtain  the  Euler-Lagrange  equations  for  (1)  by  setting  the  first  variation  to  zero  as 
/(x)  !-►  /(x)  +  Sf(x).  This  gives 

J  V(*)(/(*)  -  y{x))dx  +  X  J  tf/(i)/''"(x)ix  =  0  (2) 

for  all  variations  Sf{x)  and  hence  we  have 

xr"(*)+/(*)  =  v(*).  (3) 

The  solution  to  the  linear  equation  (3)  is  given  by 

/(»)  =  j  »■(*. »')  y{*')dx  (4) 

where  r(x,  x')  is  the  Green  function  of  (3)  obeying 

+  (S) 

Now  (5)  is  a  translation  invariant  equation  and  our  boundary  conditions  are  periodic  or  at 
infinity,  so  r(i,  x')  is  a  function  of  (x  -  x')  only.  We  write 


r(x,x')  =  r(x-x').  (6) 

For  each  different  value  of  X,  we  have  a  different  equation  (5)  and  hence  a  different 
Green's  function  which  we  denote  by  R{x,  X).  Note  from  (5)  that 


R(x,0)  =  Six).  (7) 

So  in  the  limit  as  the  confidence  in  the  data  is  complete,  the  filter  becomes  the  delta 
function.  (It  is  easy  to  see  all  this  by  taking  the  Fourier  transform  of  equation  3.) 

It  is  straightforward  to  find  the  Green’s  function  of  (5)  which  vanishes  at  plus  and 
minus  infinity.  This  is  given  by 


JKx  X)  =  — i _ e-l*l/Y^^*'*  cosT _ _ 

^  ’  ’  2X»/«  lv^X>/«  4/ 

In  Fourier  space,  the  transform  of  R{x,  X)  can  be  obtained  directly  from  (5); 


(8) 


TR(u,\)  = 


1 

1  +  Xw<’ 


(») 


so  we  can  write 

At  X  =  0  the  Green  function  goes  to  a  delta  function,  as  we  saw  in  (7). 
Define  i*  by 


Then,  for  x  >  0,  we  have 


\n*  =  1. 


(H) 
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Thus  the  extrema  of  R(x,n)  occur  at 


X  =  — nx ,  n  =  1,2... 


and  at  these  extrema,  R  takes  the  value 


So  if  li  is  small  (X  large),  the  extrema  are  at  large  x  and  correspond  to  small  values  of  R.  If 
^  is  large  (X  small),  the  nearest  extrema  occur  at  small  x  and  correspond  to  large  values  of 
R.  The  function  changes  sign  many  times. 
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Appendix  3;  The  regularization  filter  for  continuous  data  in  two  dimensions 

We  can  extend  the  results  of  the  previous  section  to  two  dimensions.  Our  generalization  of 
the  smoothing  function  /(/"{x))‘’<ix  is 


Courant  and  Hilbert  (1953)  show  that  the  Euler-Lagrange  equation  is 

=  0.  (2) 

We  write  the  regularization  of  the  two-dimensional  smoothing  problem  in  the  form 

/  / (/(x)  -  y(T))^/x  4  X  /  y (V* V/(x)V*S/(s))rfx.  (3) 

The  Euler-Lagrange  equations  of  the  combined  system  are 

W^V^V^/U)  h /(x)  =  v(x).  (4) 

This  equation  is  translation  invariant  and  so,  for  boundary  conditions  at  infinity  or  periodic, 
the  solution  can  be  written  as  a  convolution 

/(x)  =  /t-2(x)  *  »(x)  (5) 

where 

XV*V*V'‘Kj(i)4K2(i)  =  «(i)  (6) 

where  j(x)  is  the  Dirac  delta  function.  Again,  observe  that  for  small  X  the  filter  Rzix)  tends 
to  the  delta  function. 

To  solve  (6)  we -take  its  Fourier  Transform,  and  find 


T«2(«)  = 


X«*  +  f 


This  gives  a  solution 


1  f 


We  express  x  and  w  in  polar  coordinates 
So 

1  /"**  «4»wco.h(#-4) 

"*<“2;/.  i  wi-fi 

Now 


w  dw  dO. 


Jo 


where  J„  is  the  zero  order  Bessel  function.  This  gives 


n  f  \  ^  j 


that  is,  Hi  is  the  Hankel  transform  of  l/(Xw'‘  4  I).  This  integral  can  be  solved  numerically. 
The  square  Laplacian  stabilizer  leads  to  a  similar  formula,  i.e.  to 


o  /  ^  ‘  Jo[v>r)  , 

"=<*>“2/.  wTi”''"’ 


The  corresponding  filter  is  not  smooth  enough  to  allow  second  derivatives  to  be  taken,  but 
is  sufficient  for  first  derivatives.  The  integral  is  found  in  standard  texts  (for  example, 
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Gradshteyn  and  Ryzhik,  1980),  yielding 


where  kri(j-)  is  a  Thompson  function  given  by 

MX)  -  (<»l  -  fOH.)  -  fMx) .  £(-')■— --r  E 


where 


*=0  2'*"-^^+')!]' 


k=0  ^ 


and  C  =  0.5772...  is  Euler's  constant.  The  asymptotic  behavior  is 


kei(x)  i-»  "*sin(-/?x) 


where  a  and  /I  are  constants. 
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Appendix  4:  The  regularization  filter  as  an  approximation  to  the  Gaussian 

In  this  appendix  we  show  that  the  regularization  filter  approximates  the  gaussian,  both  in 
one  and  two  dimensions. 

4.1  Comparison  of  one-dimensional  filters 

We  show  that  the  one-dimensional  regularization  filter  is  an  approximate  solution  to  the 
diffusion  equation  and  therefore  approximates  a  gaussian. 

The  one-dimensional  regularization  filter  in  is  given  by 


We  expand  this  in  a  Taylor  series  in  ^  to  get 


This  expansion  is  valid  when  fix  is  small,  i.e.,  when  x  is  small  compared  to  In  this  case 
the  first  two  terms  which  we  denote  by  ft(x,n)  are  a  good  approximation  to  the  function. 
We  calculate 

(3) 

a**  2v/2  ^  ^ 


- +  0{/ia:)*, 

2^2 


which  satisfy,  to  order  the  equation 

d^ft  jdh 

dx^  ^  dft 

Thus  this  function  obeys  the  diffusion  equation, 


with  parameter  t  =  in  the  region  where  nx  is  small.  As  n  decreases,  this 

region  gets  larger  and  the  region  in  which  the  function  approximates  a  Gaussian  increases. 

This  theoretical  analysis  supports  the  numerical  results  (for  discrete  data)  which  show 
that  R  can  be  approximated  by  a  Gaussian.  Furthermore,  recalling  that  the  standard 
deviation  a  of  the  Gaussian  is  given  by  <t  =  the  analysis  shows  that  the  standard 
deviation  of  the  corresponding  Gaussian  is  X'/^. 

A  comparison  of  R  with  the  gaussian  G  =  («“*’/*"’)  can  also  be  done  directly  In  the 
Fourier  domain,  where  7R  -—  1/(1  -f  Xw^)  and  TG  =  as  shown  in  Figure  8. 

4.2  Comparison  of  two-dimensional  filters 

We  now  consider  the  two-dimensional  case.  The  regularized  filter  can  be  written  in  terms 
of  a  Fourier  integral 
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Figure  8  Comparison  of  one-dimensional  regularizing  filter  R  with  gaussian  in  Fourier 
domain:  (a)  zeroth,  (b)  first,  and  (c)  second  derivatives. 


We  perform  the  transformation  m  •-»  to  obtain,  with  ft  = 

a  f 

=  (8) 

We  expand  the  exponential  in  a  power  series 

Thus  we  have 

if)  ^  (10) 

Note  that  the  linear  term  drops  out  due  to  asymmetry  of  the  integrand.  Keeping  the  first 
two  terms  on  the  right  hand  side  of  (10),  and  denoting  this  approximation  by  we 


calculate 

•V*K  =  — 

TT 

f 

y  i-j-w® 

(11) 

and 

dh  _  1  [ 
dft  2ir  J 

1  + 

(12) 

Thus  as  before,  the  approximation  /f  j(x,X)  satisfies  the  diffusion  equation  with  t  proportional 
to  X'/'*.  The  exact  function  of  proportionality  can  be  calculated  from  (12). 


Again,  a  direct  comparison  of  the  Gaussian  with  our  regularizing  filter  Ro  is  done  easily 
in  the  Fourier  plane.  Both  filters  are  circularly  symmetric  and  therefore  depend  only  on 
the  radial  frequency  w.  A  comparison  in  Fourier  space  of  the  two  (Jimensional  gaussian 
7(7,  ^  and  the  regularizing  filter  Tll-i  l/(Xii>*  i  I)  is  shown  in  Figure  9. 


18 


Poggio 


Voorhees.  Yuille 


Figure  9  Comparison  of  two-dimensional  filter  11^  with  two-dimensional  gaussian  (7;  in 
Fourier  plane  (a)  filters  and  (b)  their  Laplacians. 
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